/*
* INSANE - Interactive Structural Analysis Environment
*
* Copyright (C) 2003-2005
* Universidade Federal de Minas Gerais
* Escola de Engenharia
* Departamento de Engenharia de Estruturas
* 
* Author's email :    insane@dees.ufmg.br
* Author's website :  http://www.dees.ufmg.br/insane
* 
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or any later version.
* 
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/
package br.ufmg.dees.insane.model.femModel.problemDriver;

import java.util.ListIterator;

import br.ufmg.dees.insane.analysisModel.AnalysisModel;
import br.ufmg.dees.insane.materialMedia.degeneration.Degeneration;
import br.ufmg.dees.insane.model.femModel.element.Element;
import br.ufmg.dees.insane.model.femModel.element.ParametricElement;
import br.ufmg.dees.insane.shape.Shape;
import br.ufmg.dees.insane.util.IMatrix;
import br.ufmg.dees.insane.util.IVector;

/** 
 * A class representing HCTProblemDriver
 *
 * @author Samir
 * @since 27/09/2006
 */

public class HCTProblemDriver extends ParametricPhysicallyNonLinearSolidMech {
	 /**
	 * 
	 */
	private static final long serialVersionUID = 1L;
	private Shape shape;
	
	/* (non-Javadoc)
	 * @see br.ufmg.dees.insane.model.femModel.problemDriver.ProblemDriver#getC(br.ufmg.dees.insane.model.femModel.element.Element)
	 */
	public IMatrix getC(Element e) {
	    int numberOfDegreesOfFreedom = e.getNumberOfDegreesOfFreedom();
		IMatrix stiffnessMatrix1 = new IMatrix(numberOfDegreesOfFreedom,numberOfDegreesOfFreedom);
	    stiffnessMatrix1.setZero();
		IMatrix stiffnessMatrix2 = new IMatrix(numberOfDegreesOfFreedom,numberOfDegreesOfFreedom);
	    stiffnessMatrix2.setZero();
		IMatrix stiffnessMatrix3 = new IMatrix(numberOfDegreesOfFreedom,numberOfDegreesOfFreedom);
	    stiffnessMatrix3.setZero();
		IMatrix G1 = new IMatrix(9,9);
	    G1.setZero();
		IMatrix G2 = new IMatrix(9,9);
	    G2.setZero();
		IMatrix G3 = new IMatrix(9,9);
	    G3.setZero();
		IMatrix finalStiffnessMatrix = new IMatrix(numberOfDegreesOfFreedom,numberOfDegreesOfFreedom);
	    finalStiffnessMatrix.setZero();
		ListIterator gns = ((ParametricElement)e).getDegenerations().listIterator();
	    AnalysisModel analysisModel = e.getAnalysisModel();
	    shape = e.getShape();    
	    IMatrix cartesianNodalCoordsMatrix = e.getCartesianNodalCoordsMatrix();
	    IVector centroid = new IVector(2);
	    centroid.setElement(0,((cartesianNodalCoordsMatrix.getElement(0,0)+cartesianNodalCoordsMatrix.getElement(2,0)+cartesianNodalCoordsMatrix.getElement(4,0)))/3);
	    centroid.setElement(1,((cartesianNodalCoordsMatrix.getElement(0,1)+cartesianNodalCoordsMatrix.getElement(2,1)+cartesianNodalCoordsMatrix.getElement(4,1)))/3);
	    
	    IMatrix cartesianNodalCoordsMatrixSubElem1 = new IMatrix(3,2);   cartesianNodalCoordsMatrixSubElem1.setZero();
	    IMatrix cartesianNodalCoordsMatrixSubElem2 = new IMatrix(3,2);   cartesianNodalCoordsMatrixSubElem2.setZero();
	    IMatrix cartesianNodalCoordsMatrixSubElem3 = new IMatrix(3,2);   cartesianNodalCoordsMatrixSubElem3.setZero();

	    cartesianNodalCoordsMatrixSubElem1.setElement(0,0,cartesianNodalCoordsMatrix.getElement(2,0));
	    cartesianNodalCoordsMatrixSubElem1.setElement(0,1,cartesianNodalCoordsMatrix.getElement(2,1));
	    cartesianNodalCoordsMatrixSubElem1.setElement(1,0,cartesianNodalCoordsMatrix.getElement(4,0));
	    cartesianNodalCoordsMatrixSubElem1.setElement(1,1,cartesianNodalCoordsMatrix.getElement(4,1));
	    cartesianNodalCoordsMatrixSubElem1.setElement(2,0,centroid.getElement(0));
	    cartesianNodalCoordsMatrixSubElem1.setElement(2,1,centroid.getElement(1));

	    cartesianNodalCoordsMatrixSubElem2.setElement(0,0,cartesianNodalCoordsMatrix.getElement(4,0));
	    cartesianNodalCoordsMatrixSubElem2.setElement(0,1,cartesianNodalCoordsMatrix.getElement(4,1));
	    cartesianNodalCoordsMatrixSubElem2.setElement(1,0,cartesianNodalCoordsMatrix.getElement(0,0));
	    cartesianNodalCoordsMatrixSubElem2.setElement(1,1,cartesianNodalCoordsMatrix.getElement(0,1));
	    cartesianNodalCoordsMatrixSubElem2.setElement(2,0,centroid.getElement(0));
	    cartesianNodalCoordsMatrixSubElem2.setElement(2,1,centroid.getElement(1));
	    
	    cartesianNodalCoordsMatrixSubElem3.setElement(0,0,cartesianNodalCoordsMatrix.getElement(0,0));
	    cartesianNodalCoordsMatrixSubElem3.setElement(0,1,cartesianNodalCoordsMatrix.getElement(0,1));
	    cartesianNodalCoordsMatrixSubElem3.setElement(1,0,cartesianNodalCoordsMatrix.getElement(2,0));
	    cartesianNodalCoordsMatrixSubElem3.setElement(1,1,cartesianNodalCoordsMatrix.getElement(2,1));
	    cartesianNodalCoordsMatrixSubElem3.setElement(2,0,centroid.getElement(0));
	    cartesianNodalCoordsMatrixSubElem3.setElement(2,1,centroid.getElement(1));

	    double[] naturalCoordsNo11 = new double[2];   naturalCoordsNo11[0]=0;         naturalCoordsNo11[1]=1;
	    double[] naturalCoordsNo21 = new double[2];   naturalCoordsNo21[0]=0;         naturalCoordsNo21[1]=0;
	    double[] naturalCoordsNo3 = new double[2];    naturalCoordsNo3[0]=1.0/3.0;    naturalCoordsNo3[1]=1.0/3.0;

	    double[] naturalCoordsNo12 = new double[2];   naturalCoordsNo12[0]=0;   naturalCoordsNo12[1]=0;
	    double[] naturalCoordsNo22 = new double[2];   naturalCoordsNo22[0]=1;   naturalCoordsNo22[1]=0;
	    
	    double[] naturalCoordsNo13 = new double[2];   naturalCoordsNo13[0]=1;   naturalCoordsNo13[1]=0;
	    double[] naturalCoordsNo23 = new double[2];   naturalCoordsNo23[0]=0;   naturalCoordsNo23[1]=1;
	    
	    IMatrix Tn1 = this.getTn(naturalCoordsNo11,naturalCoordsNo21,naturalCoordsNo3,cartesianNodalCoordsMatrix,0);	    
	    IMatrix Tn2 = this.getTn(naturalCoordsNo12,naturalCoordsNo22,naturalCoordsNo3,cartesianNodalCoordsMatrix,12);	    
	    IMatrix Tn3 = this.getTn(naturalCoordsNo13,naturalCoordsNo23,naturalCoordsNo3,cartesianNodalCoordsMatrix,24);	    
	    while(gns.hasNext()) {
		    Degeneration degeneration = (Degeneration) gns.next();
	        double[] naturalCoords = degeneration.getRepresentation().getNaturalcoords().toDouble();
	        IMatrix matrixC = degeneration.mountC();
		    IMatrix G = new IMatrix(9,9);   G.setZero();
		    IMatrix auxG = new IMatrix(9,3);   auxG.setZero();
		    auxG.mulTransposeLeft(getFi(naturalCoords),matrixC);
		    G.mul(auxG,getFi(naturalCoords));
	        
	        // subElement1
	        IMatrix derivedShapeFunction1 = shape.getDerivedShapeFunction(naturalCoords, cartesianNodalCoordsMatrixSubElem1);
			double weight = degeneration.getRepresentation().getWeight();
			double jac1 = analysisModel.getJacobianTransformation(derivedShapeFunction1,cartesianNodalCoordsMatrixSubElem1);
			double coeff = e.getElmCoefficient();
			IMatrix matrixG1 = new IMatrix(9,9);   matrixG1.setZero();
			matrixG1.add(G);
			matrixG1.setScale(coeff*weight*jac1);
			G1.add(matrixG1);
			
			// subElement2
	        IMatrix derivedShapeFunction2 = shape.getDerivedShapeFunction(naturalCoords, cartesianNodalCoordsMatrixSubElem2);
			double jac2 = analysisModel.getJacobianTransformation(derivedShapeFunction2,cartesianNodalCoordsMatrixSubElem2);
			IMatrix matrixG2 = new IMatrix(9,9);   matrixG2.setZero();
			matrixG2.add(G);
			matrixG2.setScale(coeff*weight*jac2);
			G2.add(matrixG2);

			// subElement3
	        IMatrix derivedShapeFunction3 = shape.getDerivedShapeFunction(naturalCoords, cartesianNodalCoordsMatrixSubElem3);
			double jac3 = analysisModel.getJacobianTransformation(derivedShapeFunction3,cartesianNodalCoordsMatrixSubElem3);
			IMatrix matrixG3 = new IMatrix(9,9);   matrixG3.setZero();
			matrixG3.add(G);
			matrixG3.setScale(coeff*weight*jac3);
			G3.add(matrixG3);
	    }
	    IMatrix auxSM1 = new IMatrix(Tn1.getNumCol(),G1.getNumRow());   auxSM1.setZero();
	    auxSM1.mulTransposeLeft(Tn1,G1);
		stiffnessMatrix1.mul(auxSM1,Tn1);

	    IMatrix auxSM2 = new IMatrix(Tn2.getNumCol(),G2.getNumRow());   auxSM2.setZero();
	    auxSM2.mulTransposeLeft(Tn2,G2);
		stiffnessMatrix2.mul(auxSM2,Tn2);

	    IMatrix auxSM3 = new IMatrix(Tn3.getNumCol(),G3.getNumRow());   auxSM3.setZero();
	    auxSM3.mulTransposeLeft(Tn3,G3);
		stiffnessMatrix3.mul(auxSM3,Tn3);

		finalStiffnessMatrix.add(stiffnessMatrix1);
		finalStiffnessMatrix.add(stiffnessMatrix2);
		finalStiffnessMatrix.add(stiffnessMatrix3);
		
		return(finalStiffnessMatrix);
	}
	
	private IMatrix getFi(double[] naturalCoords){
		double e1 = naturalCoords[0];
		double e2 = naturalCoords[1];
		double e3 = 1-e1-e2;
		
		IMatrix fi = new IMatrix(3,9);
		fi.setZero();

		fi.setElement(0,0,e1);   fi.setElement(0,1,e2);   fi.setElement(0,2,e3);
		fi.setElement(1,3,e1);   fi.setElement(1,4,e2);   fi.setElement(1,5,e3);
		fi.setElement(2,6,e1);   fi.setElement(2,7,e2);   fi.setElement(2,8,e3);
		
		return fi;
	}
	
	private IMatrix getTn(double[] nC1, double[] nC2, double[] nC3, IMatrix cN, int k){
        IMatrix secondDerivedShapeFunction1 = shape.getSecondDerivedShapeFunction(nC1, cN);
        IMatrix secondDerivedShapeFunction2 = shape.getSecondDerivedShapeFunction(nC2, cN);
        IMatrix secondDerivedShapeFunction3 = shape.getSecondDerivedShapeFunction(nC3, cN);
	    IMatrix Tn = new IMatrix(9,12);   Tn.setZero();
        for (int i = 0; i < 12; i++) {
			Tn.setElement(0,i,secondDerivedShapeFunction1.getElement(0,(i+k)));
			Tn.setElement(3,i,secondDerivedShapeFunction1.getElement(1,(i+k)));
			Tn.setElement(6,i,(2*secondDerivedShapeFunction1.getElement(2,(i+k))));

			Tn.setElement(1,i,secondDerivedShapeFunction2.getElement(0,(i+k)));
			Tn.setElement(4,i,secondDerivedShapeFunction2.getElement(1,(i+k)));
			Tn.setElement(7,i,(2*secondDerivedShapeFunction2.getElement(2,(i+k))));

			Tn.setElement(2,i,secondDerivedShapeFunction3.getElement(0,(i+k)));
			Tn.setElement(5,i,secondDerivedShapeFunction3.getElement(1,(i+k)));
			Tn.setElement(8,i,(2*secondDerivedShapeFunction3.getElement(2,(i+k))));
		}
        return Tn;
	}
}
